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I. INTRODUCTION 


Consider a nonrigid structure submerged in a fluid. 
Pressure in the fluid may induce response motions of the 
structure. Conversely, motion of the structure affects 
the pressure field. Such a system is called a Coupied 
elasto—hydrodynamic system. 

Many practical engineering problems belong to this 
fereeoey A reservoir-dam system subjected to earthquake 
motion, oscillations in rocket fuel systems, a vibrating 
ship's hull, and structures subjected to blast loadings due 
to explosions are familiar elasto-hydrodynamic problems. 

The particular elasto-—hydrodynamic system under consid- 
eration in this paper is a submarine subjected to a blast 
wave resulting from an underwater explosion. The submarine 
is modeled as a ring-reinforced, neutrally buoyant, circular 
eylinder with rigid bulkheads. It is assumed that the blast 
wave may be modeled as a ane acoustic wave with wave-front 
parallel to the longitudinal axis of the submarine. 

The interaction of an acoustic wave with an obstacle has 
been the subject of many investigations. Sette [1] first 
computed the pressure distribution on a rigid cylindrical 
surface due to a step-pulse. The response of a cylindrical 
shell to a plane step (or decaying) shock wave has been 
treated by various authors [2-5]. 

Due to the complex nature of the coupled response a 


realistic analytic formulation of this problem is intractable. 
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However, numerical methods, particularly the finite element 
method, provide a satisfactory means of solution. 

The finite element formulation of structure-fluid 
dynamics was first introduced by Zienkiewicz, Irons ang 
Nath [6]. They demonstrated that the coupled natural 
frequencies and mode shapes of systems could be accurately 
determined by this method. 

It is the purpose of this paper to demonstrate that the 
finite element method, as forwarded by Zienkiewicz and 
Newton [7], is applicable to the solution of the interaction 
of .a blast wave with a rigid structure. Using the finite 
element formulation, a superposition theorem is developed 
which affords a method of applying the rigid cylinder results 
to a realistic, but simple, three-dimensional model of t 


submarine. 
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It. THE FINITE ELEMENT METHOD OF SOLUTION 


A. THE STRUCTURE DISCRETIZATION 

The finite element discretization of the structure 
problem is well known and is described in Zienkiewicz's 
text. [8]. The discretization process results in an 


assembled system of equations of the form 
[M]{6} + [c]{S} + [K]{6} = {R} . (1) 


In the above equation [M], [C] and [K] are respectively the 
mass, Gamping, and stiffness matrices of the structure 
calculated in the proper, consistent, manner, {6} is a vector 
of nodal displacements, the dots indicate time differentiation, 
and {R} is a vector of generalized nodal loads. 

It is convenient to divide the generalized nodal forces 


into two parts (after Zienkiewicz and Newton [7]) such that 
tro LE PD (2) 

The vector {F} is due to external forces and {P} is due to 

the fluid pressure on the interface. For the ge structural 

node 

t= 4 N's as 5 (3) 

S 
where N,! is the appropriate shape function defining the 


displacement pattern in the direction normal to the boundary, 


a3 
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p is the pressure on the structure-fluid interface S, and 


the integration is over the interface. 


B.. THE FLUID DISCRETIZATION 
From first principles it can be shown [6,9] that the 


wave equation 
sles 
ig eee ie ae Ce, 
c 


= is the Laplace operator, p is the hydrodynamic 


where V 
pressure in excess of static pressure and c the acoustic 
velocity of propagation, together with the necessary 
boundary and initial conditions are the governing relations 
for the hydrodynamic response of the fluid. 

The spatial discretization of the fluid problem is 
accomplished by a direct application of the Galerkin weighted 


residual process. If at any instant of time p is approximated 


by 
Dp =a. G. (summaeion on J = l....,m). , CS) 


th 


Where m is the number of nodes, N. the j shape function 


J 
chosen to describe the spatial variation of p over the region 


and P, is a set of nodal pressure values which are time 


th 


dependent, then the i weighted residual equation is 


2 eral se ie 
4 N LV (NsP 5) = ie (N5B,)] dR Pee (6) 
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where N, is the shape (weighting) function [8, Ch. 3] of 
the Galerkin process and R is the fluid region under 
consideration. Employing Green's theorem, the weighted 


residual equations may be transformed to 
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This system of equations can be written in matrix form as 


cigs LEY ey (8) 


The matrices [Q] and [H] are the assembled "inertia" and 
"stiffness" matrices of the fluid. They are defined on the 


element level by: 


1 
q.; =—5 JS N,N, aR (9) 
1J ef R a a) ae? 
e 
oN. ON. oN. oN oN. oN 
= eeteeey Seema pe A. J 
Ds j [ x Ox dy oy 0Z OZ i oe (10) 
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where R denotes the element region. 


The vector {B} results from the surface integral in (7) 


and can be written for the 4th surface node as 
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where n is a local coordinate in the direction of the 
outward normal. By means of this vector it is possible to 
apply the boundary condition on the normal pressure gradient, 


Bee a5 


ree oe Py ey 


where Vn is the fluid particle velocity normal to the 
boundary and p is the fluid density. It is noted that for 
a fixed boundary ue = 0 and the corresponding bs =0. 3 

Ay nonreflecting' or radiation boundary condition, 
necessary to model numerically an infinite fluid region, was 
developed by Zienkiewicz and Newton [7]. They demonstrated 
that for a plane wave normally incident upon a plane boundary 


the condition to be satisfied by p on the boundary S is 


@ 


20s ag hee 
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The corresponding boundary integral (11) is 
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This leads to an additional "damping" term [D]{p} in the 


Matrix equation (8), i.e., 


(Ql{p} + [D]{p} + CH]{p} = {B} . (15) 
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{[D] is the assembled fluid "damping" matrix defined on the 


element level as 


goes (16) 


=i 
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where Se denotes the external boundary of the element. 


C. THE COUPLED STRUCTURE~FLUID .PROBLEM 

The vectors {P} of the force term of equation (2) and 
{B} of equation (15) determine the coupling of the fluid 
and structure in the discretized problem. Recall from 
equation (3) that the nodal forces on the structure due to 


the pressure are 
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where [L] is the assembled matrix of element contributions 


toe = 2 N,N, ds. . (18) 


oe Sis oS aa! 8 ; (19) 


a7 


Pet hart ( Hite) 8 Shi tt 


hears igen 


A al Be (ae: eo 5 ay ong 














(ax) ee neta 
«Yeattele eff Io Yxsbovod Esivtadxe: aa t) 


watEors casein 

bee (2) not'Seupe to artey sSxot edt to esp" 
‘ESE dad ta Antiques ad? scettrresob (et) § 
nor? flasen’ “yma ety bas ttonserb ond wee 

ae we yeeeewere S439 no Bsoto?” tehon’ aay Pe 


ee a 


=a 


ya Tee We) = a Ge : tay) ay 


Hl oy y 
i $2 “ 


: 7 - 


« (at ni) =| 


(TI) 


i ibe 
? ' . , i ee 
knotzudiyviwod Jremuls ‘lo xPdaw belduanen ste ae ‘tay 


(Bi) . 8b wn YOu? 
) ‘ t ; ts? 
= 


if 
(SY aontaape nate | 


2 'yva- * ee “a 








The forcing term {B} in equation (15) becomes 


B,\ 2 /-PF NIC N.' 68 
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or 


(5 kese=alS HéiLte, (20) 


with [S] the assembled matrix of element contributions 


S45 = 0 a NAN, mae tt (21) 
Observe that 
aa 


The complete elasto-hydrodynamic system has been 
formulated and can be summarized. Rewriting equation (1) 
wich (2) and (17), 


[mM]{S} + (c]{é} + [K]{6} = [L]{p} + {F} (23) 


describes the structure behavior. Similarly, equation (15) 


with (20) and (22) can be written as 
bodip se (bite) + Liliols= =p. [ul te}. (24 


This governs the fluid response. 
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The complete system of equations can be written as 


Bie we ST: ee 


Using general finite element techniques it is possible 
to generate and subsequently solve this unsymmetrical system 
of equations for the submarine-fluid problem under considera- 
tion. However, an alternative method based on this general 
formulation will be utilized to effect a solution using the 


following superposition theorem. 


D. SUPERPOSITION THEOREM 

It will be shown that the coupled elasto-hydrodynamic 
problem formulated above may be solved in two stages and that 
the structure response thus found is the same as in the 
complete unseparated problem. In the first stage the structure 
is constrained against responding to the loading imposed by 
the incident pressure pulse. In the second stage the structure 
and the fluid are initially at rest. External loads are 
applied to the structure to simulate the loading by the 
pressure pulse of the first stage and the structure and fluid 
responses are found. 

1. The Complete Problem 

Equations (23) and (24) govern the structure and 


fluid response subject to the initial conditions 
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{8(G) he Sener (26) 


{p(0)} fait wes feCOly => BT vs (273 


where {a} and {b} are given vectors. It is assumed that 
a ae a ae (28) 


The solution of the complete problem requires that equations 
(23) and (24) with equation (28) be solved for {6(t)} and 
{p(t)}, subject to the initial conditions, equations (26) 
and (27). 
2. wsvape One 
* : (aly) \ ~ re 1 
The dependent variables will be {p‘*’(t)} and {6° ’(t)}. 
Recall that in the first stage the structure is constrained 


against motion, i.e., 
Heer) Pye ok. (29) 
For the fluid, equation (24) becomes 
Eadge 73 + foie} + tayie y= 0. , (30) 
with the initial conditions 


nd tale fe (0d = tek (31) 
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Equation (23), governing the structure response, will be 


satisfied by (29) if it is required that 

eS Sete ty (32) 

For the first stage one need only solve equation (30) 
for ip‘?) (t)} subject to the initial conditions of equation 
ied ug cp), is subsequently found from equation (32). 

3. Stage Two 

The system is initially at rest and the structure is 
externally loaded to simulate loading by the pressure pulse. 
The dependent variables for this stage are p62) (4)} and 
ca°?) (4). For the fluid 
Iotins?\) + Epitps?’) 4 pee) = -pcoitis + , 43a) 
with the initial conditions 

tp6?2(0)} = (822 (0)} =o. (34) 
For the structure 
cmits’??) + po7(8(?s + ceqge’??> = crim '?)) + trf?)), (35) 
with the initial conditions 


Ge ny) ee =) (01 = 0 (36) 
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Specify further that 
(HE) ment gon . (37) 


For the second stage of the solution it is necessary 
to solve equations (33) and (35) for tp?) (¢)} and (62) (¢)} 
with equation (37) subject to the initial conditions of 
equations (34) and (36). 

4, Superposition 


It is now asserted that 


rice = es ee eee (38) 
and 


eaqene a te (eh eA ey (39) 


In virtue of the conditions imposed on the first 
stage pressures {p(y} and displacements (661) (453 : 
equations (29) - (32), and on the second stage pressures 
(p62) (¢)} and displacements (6 (2) (¢)} , equations (332) = (37); 
it follows that {p(t)} and {6(t)} as given by equations (38) 
ana (39) satisfy the governing equations for the complete 
_ problem, equations (25) - (28). This establishes the super- 
position theorem for the coupled elasto-hydrodynamic problem. 
Implementation of the superposition theorem does not 
require that the rigid body pressure distribution or the 


structure response be obtained by the finite element method. 
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It is noted that for the submarine structural model under 
consideration a three-dimensional analysis is necessary. 
However, the superposition theorem permits the rigid body 
pressure distribution to be obtained from a two-dimensional 


analysis. 


23 


Vitae 


ad he 


Li a, 





IIIT. WAVE PROPAGATION BY THE FINITE ELEMENT METHOD 


To evaluate numerically the interaction of a wave with 
a structure it is necessary to propagate the blast wave 
through the contiguous fluid region. Before proceeding with 
the complicated two-dimensional interaction of a blast wave 
with a rigid cylindrical structure, the characteristics of 
the one-dimensional wave propagation are studied. The 
techniques developed and utilized for the one-dimensional 


case are subsequently employed in the two-dimensional study. 


A. ONE-DIMENSIONAL WAVE PROPAGATION 

The numerical results obtained by the finite element 
method for the one-dimensional wave propagation problem are 
easily compared to well known theoretical values and thus 
provide insight into the behavioral characteristics of the 
techniques employed. 

The physical situation to be considered is a semi-infinite 
fluid strip of unit cross-sectional area. Through this region 
a pressure wave of known shape is propagating with acoustic 
velocity c. A rigid structure is located at x = 0 and the 
radiation boundary condition is applied at x = L. Equation 


(24), with {6(t)} = 0, 
[a]{p} + (D]{b} + [H]{p} = 0 , (40) 


governs the fluid response for the case. The initial 
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conditions {p(0)} ana {p(0)} are determined by the wave form 
under consideration. 
1. Blast Wave Model 

A blast wave, which is generated by an underwater 
explosion, is characterized by an extremely rapid rise in 
pressure to a maximum value p*, followed by decreasing 
pressure. For the purpose of investigating blast wave 
propagation it would appear that a step pulse of finite 
duration would be a satisfactory wave model. However, it is 
evident that to model a moving discontinuous function by a 
finite fixed spatial discretization is not possible. Conse- 
quently, for the purpose of this investigation the blast wave 
is modeled as a ramp rise to a constant value, i.e., a ramp- 
step wave. Although the ramo-step wave is of primary interest, 
this investigation is not limited exclusively to that wave 
form and illustrative results of tests made using different 
wave forms are also given later. 

2. Element Selection 

To study a significant portion of a blast wave, which 
is propagating at the acoustic velocity, requires a spatial 
domain of considerable extent. This is necessary to avoid 
spurious effects caused by reflections from the (artificial) 
boundaries. Accordingly, the required linear dimensions of 
the region are approximately proportional to the time interval 
during which the structure-fluid response is sought. It is 
evident that, if a fine subdivision of the one-dimensional 


fluid strip is required to represent the wave form under 
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consideration, the number of degrees of freedom in PERN 
(40) may become prohibitively large when the same subdivision 
is applied to the two-dimensional region. 

Two types of finite elements, linear and cubic, are 
considered for the spatial discretization of equation (40). 
As their names imply they respectively provide for a linear 
and cubic variation of the pressure within an element. When 
propagating smooth wave forms represented by a significant 
number of elements, e.g., p(x) = p* sin? (1 x/L) or 
p(x) = p*cos (1tx/L), both the linear and cWhie ~okemeues 
produce excellent numerical results. However, in choosing 
the type of finite element best suited to represent spatially 
a propagating blast wave it is necessary to consider: a) the 
general shape of the wave; and b) two instantaneous finite 
element representations of the wave front as it propagates 
across a region. 

As established ees the general shape of the wave 
of primary interest is the ramp-step wave. For the second 
consideration, assume that the node spacing Sn is yequadl /ror 
both element types. Due to the character of the linear and 
cubic elements note that three linear elements are required 
for each cubic element in the spatial representation of the 
wave front. Now consider a ramp wave front that at one 
instant in time is represented spatially by four nodal values, 
i.e., the ramp rise spans one cubic element. At this 
particular instant both element types are capable of an exact 


representation of the ramp wave front. For the second instant 
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in time, let the wave advance a distance equal to one node 
space. At this instant of time the linear elements continue 
to provide an exact representation of the wave front. 
However, when the wave is in this spatial position relative 
to the cubic elements, they are required to fit a cubic 
polynomial through a function which has a discontinuous first 
derivative. This results in a distortion of the wave front. 
The above qualitative discussion explains the 
observed superior performance of the linear element when 
used in the spatial discretization of equation (40). It also 
serves as a reminder that higher order finite elements per 
se do not guarantee superior results and may indeed lead to 


poorer results. 





3. Wime Integration 
The time integration of the governing equation (40) 
is computed using two finite differencing techniques, the 
Newmark-8 method as adapted by Chan, Cox and Benfield [10] 
and the Houbolt method [11]. Both provide a numerical method 


of determining the time dependent response of a differential 


equation of the form 
mw + aw t+ kw = f(t) : (41) 


where m, dad and k are constants. The difference equation 
proposed by Chan, et al. [10] is 
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where h is the time step and 8 is a parameter which may have 
any value from 0 to 1/4. The value B = 1/12 is chosen for 
calculations using the stability requirements established by 
Newmark [12]. 


The starting equation is 


3 
(m+ Ba + phok)w, = (m+Ba- - py - - 8) PH wv, 


Die 
+ (m - (F -e 2S) hig + gh°r, 


+( 5-8) + G- 8) Syne, (43) 


where Wo and Wo are initial values. 
The Houbolt method [11] is developed from a consider- 
ation of a cubic curve that passes through four successive 


ordinates to obtain the difference equation 


(2m + ae 


2 iJ 3 
dhst kh Ds = (5m + 3dh )w - (4m + 5 dh)w. 
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The starting procedure for the Houbolt method employed 
in the one-dimensional propagation problem is as follows: 
assume that the wave form is known at t = 0 and that the wave 


is propagating undeformed through the fluid region with the 
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acoustic velocity c. The pressure distribution p can be 


described as 
p(x,t) = f(x + ct) (45) 


for a wave moving in the negative x-direction. It is apparent 
that once the wave form is specified at t = 0 as a function 
of x the required preceding values of pressure may be calcu- 
lated directly from equation (45) and thus provide exact 
starting values for equation (44). 

It is noted that the Houbolt method Me Gheore seas 
stable for all step sizes, while the Newmark-8 method becomes 
unstable if the limiting step size established by Newmark [12] 
ue exceeded. 

_As an example of the results obtained by the methods 
outlined above consider a fluid strip in which the acoustic 
velocity c = 5000 ft/sec. The region of interest extends 
HeuO fc. from/a-rigid structure located iat x°= 0 and is 
divided by 31 equally spaced nodes. The integration time 
step h = 0.01 sec. At time t = 0 the wave form has the shape 
and position in the region shown in Fig. 1. As the wave 
strikes the rigid structure at x = 0 it is reflected and 
propagates to the right where the radiation boundary condition 
is applied. 

After 100 time steps (t = 1.0 sec.) the wave will be 
positioned as shown in Fig. 2. The corresponding numerical 
results are shown for the Newmark-8 integration with cubic 


elements and the Houbolt integration - linear element 
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Fig. 1. Pressure-vs.-distance at t = 0 
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combination. The other combinations of integration techniques 
and elements give essentially identical results. 

It is evident that to determine the optimum integra- 
tion technique it is necessary to establish criteria to judge 
the results of tests. The criteria chosen are: a) the 
pressure-time history at the structure-fluid interface 
(impulse) and b) retention of wave shape after propagation 
through the fluid. 

If the rise time tg) of the ramp is chosen as the 
parameter necessary to describe the wave front, criteria for 
selecting the time h and the node spacing S, can be developed. 

Assume that to describe the slope of the ramp rise 
at least two linear elements are required, i.e., 


2s 
t > (46) 


e Cc 
Intuitively this restriction on the node spacing seems 
reasonable because it will provide that at least two spatial 
locations in the fluid strip are used to describe the slope 
of the ramp. 
It also seems reasonable to restrict the time step 
in such a way that the wave front cannot cross an element in 


one-step, i.e., 


Ss 


n ! 
Lea 7 (47) 
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The above restrictions are proposed as a general 
guide in the initial selection of the nodal spacing and the 
integration time step. They are not intended to establish 
restrictions which cannot be violated. As shown below, 
Satisfactory results may be obtained when the restriction on 
t (46) is violated. 

Using the criteria established above to judge the 
effectiveness of the integration technique, numerous tests 
were made to determine which time integration technique is 
better Bh ecd<fon the purpose of propagating a ramp-step wave. 

Consider, for example, a triangular shape wave that 


has the following characteristics: 


s 
= een, 
Be etl Cc > 
Sh 
oc ane é 
Sh 
te aay 5 


where tT, is the time required for the wave to decay from the 


a 
maximum value p* to zero. The pressure-time histories at the 
rigid structure-fluid interface obtained by both integration 
methods are shown in Fig. 3. The result shown in Fig. ae 
obtained by employing a cubic finite element discretization 


with the Houbolt integration technique for the same wave form. 


Consider a ramp-step wave with the characteristics: 
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Pressure-vs.-Time at the Rigid Structure- 
Fluid Interface (ramp-step wave) 
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The resulting pressure-time histories at the Sig eee eee 
interface are shown in Fig. 5. From this and other results 
it is concluded that the Houbolt integration technique is 
better suited for the propagation of the blast wave through 


a fluid region. 


B. TWO-DIMENSIONAL WAVE PROPAGATION 

Based upon the results of the one-dimensional study, the 
linear finite element and the Houbolt integration technique 
are selected to evaluate numerically the pressure-time 
history at the rigid cylindrical structure-fluid interface 


resulting from a ramp-step wave. The geometry of the problem 


= shown in Fig. 6. 





Fig. 6. Two-dimensional geometry 
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The rigid cylindrical structure, of radius a, is centered 
at the origin of the x,y coordinate system. The diagram- 
matically shown ramp-step wave is propagating through the 
fluid in the negative x-direction with the acoustic velocity 
ec. Due to symmetry only the upper half plane is shown and 
the rigid body boundary condition is applied along the 
x-axis. The additional necessary boundary and initial 
conditions are specified below. It is convenient to define 
a nondimensional time parameter t = ct/a, where t = 0 is the 
time of incidence of the wave at r = a, 6 = 0 

1. Boundary Conditions 

The radiation boundary condition as developed earlier 
is valid only for a plane wave normally incident upon a 
Plane boundary. Clearly the wave reflected from the rigid 
cylinder will not be normally incident upon the far fluid 
boundaries. It is necessary that these boundaries be 
located a sufficiently large distance from the structure to 
insure that the results at the interface are not affected. 
This investigation and Sette's [1] results both indicate that 
the transient response of the blast wave-cylindric structure 
Bees cin is essentially completed in 8a/c time units after 
incidences at r = a, 6 = 0% As a consequence the far fluid 
boundaries are located sufficiently far from the structure 
to prevent any boundary reflected waves from reaching the 
structure in the time specified. For programming convenience 
the rigid body boundary condition is applied to the boundary 
at y = L and the radiation boundary condition is applied at 


x = +L. 
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A significant characteristic of the radiation boundary 
condition was observed during the two-dimensional study. It 
was observed that if the constant value portion of the ramp- 
step wave is applied to the boundary nodes at x = L for three 
successive time increments, then the radiation boundary 
condition continued to let the wave enter the fluid region of 
interest. In addition the interaction of the wave reflected 
from the structure with the boundary does not affect the 
incoming wave. This observation permitted a significant 
reduction of the fluid region under consideration, resulting 
in a proportionate reduction in the number of degrees-of- 
freedom required for the spatial discretization of equation (40). 

2. Two-Dimensional Mesh Consideration 

In considering the two-dimensional spatial representa- 
tion (mesh) of the fluid field by the finite element method, 
it is necessary to incorporate both the physical geometry and 
the wave propagation aspects of the problem. As described in 
the one-dimensional study the integration time step is 
logically limited by the node spacing. Equation (47) governs 
this dependence. The physical dimension of the cylindric 
structure, when considered in the context of the advancing 
wave front, also restricts the nodal spacing. Clearly. at 
the node spacing s,, and structure radius "a" are chosen to 
be of the same order of magnitude and the upper limit on the 
time step is utilized, the wave front could advance across 
the structure ina few time steps and the pressure-time 


history at the interface would be of suspect value. It is 
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therefore necessary to produce a mesh which affords a degree 
of fineness near the structure compatible with the integration 
time step and wave characteristics. Additionally, computer 
storage and processing time restrictions must be considered. 
After extensive experimentation the mesh shown in 
Fig. 7 was developed. Due to geometric symmetry about the 
y-axis only the region of the first quadrant is shown. This 
mesh has the following notable characteristics: 

a) Elements which form the cylindrical structure boundary 
are ‘'special' 5 noded, linear-parabolic elements. At 
the interface the 'special' elements represent the 
rigid boundary by a parabolic are and provide for a 
transition to the linear elements in the fluid field. 
The interface nodes are equally spaced at angular 
invervals: of 7.5°. 

b) The 'near fluid' region is characterized by a radial 
node spacing equal to the 'far fluid' node spacing and 
an angular node spacing of 15°. 

ec) The transition from the cylindrical 'near fluid' region 
to the rectangular 'far fluid' region is accomplished 
with a minimum number of distorted and degenerate 
elements. 

d) The uniform rectangular 'far fluid' region is charac- 
terized by node spacing 0 5 a. The appropriate 
boundary conditions as discussed above are applied at 


x = +6a and at y = 6a. 
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The overall characteristic of the mesh is element 
size uniformity. Results indicate that such uniformity in 
both spatial directions is desirable for satisfactory two- 
dimensional wave propagation. 

3. Integration Starting Procedure 

The Houbolt integration starting procedure, equation 

(45), used for the one-dimensional case assumed a general 


solution of the wave equation (4) of the form 
p(xe,t) = P(x Fev + g(x = ct) (48) 


with g(x-ct) = 0. This solution is valid assuming that 
f'(x) and f"(x) exist. However, the ramp-step wave has 

a discontinuous first derivative and singular higher 
derivatives. To avoid this inconsistency assume that the 


ramp wave front takes the form of a cubic curve at t = 0. 


% 
p= 2+ ox! + a(x")? , (49) 
with = = q + 3B(x")* (50) 
Where x' = x - Xo 
= dp oo ss 1 . 8 
If p = p* and aes O.at.x'! = €.as shown in. Fig. 


it is then possible to solve equations (49) and (50) for the 


unknown coefficients 
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Fig. 8. Cubic Wave Front 


= 5p" pe 
a = : and BF-—= , 
TE 4e3 
or 
pedtpeci+ 3% -3 (7) . (52) 


By substituting into (48) and differentiating with respect 


to time 


p - 2B oa-(%)*) . (53) 
t=0 
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Equations (52) and (53) determine the initial conditions 
necessary for a proper formulation of the boundary value 
problem in the region of the cubic ramp wave front. 

To incorporate the initial conditions into a starting 
procedure for the Houbolt integration technique consider a 


Taylor series expansion about t = 0 
Pecapg tpgtriDgt ., (54) 


where Po is p : 
t=0 


If h is the time step 


e he 
Poa Pp pt. a Po > 
(55) 
, ne“ 
Bee Goes aE On: 
The vector {py} is deduced from equation (40), i.e., 
{pp} = Cal (-[c]{By} - [H]{pg}) . (56) 


Thus, equations (55) with equations (56), (52) and (53) 
provide a consistent starting procedure for the two- 
dimensional propagation problem employing the Houbolt 
integration technique. 

To minimize errors in the numerical solution of the 
pressure-time history, the wave front is positioned in the 
fivid so that at t = 0 the cubic-ramp has its initial (leftmost) 


peedewaa P= a, 6 = 0°. 
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Consider the following integration, wave and mesh 


characteristics: 
hos gees 
Th = een : 
an = awe 


Implementation of boundary conditions, mesh and 
starting procedures described above results in the pressure- 
time histories at 6 = 0° and 6 = 180°, shown in Figs. 9 and 
10. The solid line represents Sette's [1] results displaced 
along the time axis to coincide with the incidence of the 
constant portion of the cubic ramp-step wave. Additional 
curves are given for 6 = 30°, 60°, 90°, 120°, 150° in 
Appendix A. No comparison is available for the curves given 
for 120° and 150°. Tabulated values of pressure vs. time at 
each structural node are given in Appendix B. Figures 11 and 


12 provide a comparison of results when the time step is 


OuHO A 


increased to h = a2 and the rise time decreased to te = 2 


3¢ 
The results shown in Figs. 9 and 10 and Appendix A are 
in good agreement with those calculated by Sette ie Ne em oe 
agreement establishes that the Houbolt [11] integration 
pees and the linear finite element are a suitable 


combination for evaluating numerically the interaction of a 


blast wave with a rigid cylinder. 
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Fis. 9. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (6 = 0°) 
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Fig. 10. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (6 = 180°) 
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Fig. 12. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (6 = 180°) 
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IV. THE SUBMARINE STRUCTURAL MODEL 


A simple three ta seco ott rockon submarine structural 
model is developed. This model demonstrates the implementa- 
tion of the superposition theorem using the pressure-time 
history obtained from the two-dimensional finite element 
solution of the interaction of a blast wave with a rigid 
cylinder. 

A single, ring reinforced section of a submarine, 
terminated at both ends by bulkheads, is modeled. The effects 
of the stiffening rings are included by treating the ey 
as orthotropic. Each bulkhead is rigid in its own plane. / 
The added mass effect due to the surrounding water is 
included by considering the water to be incompressible and 
using strip theory. 

The three modes considered are: 1) rigid body mode, 

2) shear beam mode and 3) buckling mode. The equations of 


motion take the form 
[M]iq} + [K]{q} = {f} , (57) 


where [M] and [K] are the structural mass and stiffness 
matrices respectively, {aq} the vector of generalized displace- 
ments ana {f} the vector of generalized forces. 
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Fig. 13. Submarine Structure Geometry 


The following symbols are used to describe the structure 
model: 


= shell radius, 
= shell thickness 
= bulkhead spacing, 


n 


= radial displacement 
tangential displacement 

= axial displacement 

= Young's modulus of elasticity, 
= shear modulus of elasticity 


<<. Gi Pies <j Ss 2siny © 
i 


= Poisson's ratio. 
Values of the structural parameters used in calculations 


are specified in Appendix C. 
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The generalized modal displacements are defined as: 


1) Rigid body mode, 
u = -q, cos Og eS qy gin 8, w=0. (58) 


2) Shear beam mode, 


aie sin > eos 6°, v= qd, sin —- sin 8 . 
wee 0 4 (59) 
3) Buckling mode, 
Pe ivAa oe 1 : 
uU= Ton sin ZT cos 28 : at Snes q3 sin > sin.26. , 
w= A d3 sin oye eos 28) , (60) 


The parameter X is evaluated by minimizing the modal strain 
energy. 

In the rigid body mode and the shear beam mode, each 
cross-section (z = constant) moves as a rigid body in the 
athwartship direction. In the buckling mode a cross-section 


is deformed, but the tangential strain Eg = O. 


A. ADDED MASS EFFECTS 

Conceptually the effects of fluid pressure Pa resulting 
from the radial motion of the structure should be assigned 
to the generalized force vector. It will be shown that the 
corresponding components are linear functions of qdy> do and 
d3 so that it is appropriate to include these effects by 


augmentation of the mass matrix. 
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Consider first a shell segment of infinitesmal length 


(dz) for which the radial displacement is 


u =u, cos 6 + uy cos 26 , (61) 
= NZ 
where U, = -d) - dy Sin 7 e (62) 
and us, = -qg sin® ae (63) 
2 3 vy " 


The pressure Pin exerted on the hull as a result of 
radial motion is determined using potential theory. For an 


incompressible fluid 


¥ ba. =. 0s (64) 
Let P= By Pa ‘ (65) 
1218 Py = Po = cos 6 f(t) : (66) 
ac 
and P5 = Pg a COS 26 g(t) 5 (67) 
r 


then Pe satisfies equation (64). 


‘From Euler's equation of motion for an inviscid fluid, 


ap a 
a - pu, cos ae (68) 
Py = pau, cos Bin ty (69) 
r= 
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where op is the fluid density. Similarly, 
Pa| = oe Uy cos 26 ; (70) 
Therefore, the pressure is given by 


pa(u, cos 6 + - u, cos 26) 


P| = [Pp] + Polyea 2 


r=a 
(71) 


Let ow" be the work done by the pressure axe during a 
virtual displacement 6éu. 
Then 
So = eae p. Su ade) (72) 
Substituting 


éu = duy ecs 6 + dus GUSV 28 | | is 


and performing the indicated integration gives 


éW,! = ~mea® (uy éu, + = Us 6Uy) dz. (72) 
If the substitutions 

du, = -64, - 645 sin > r (73) 
and 

Su, = ~644 sin sae : (74) 
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are made and the integration over the length carried out 


we have 
om if 
m =e 0 qy 
if 2m m a 
3m 


where m = oma’? (the displaced mass). The added mass matrix 


is 


om 
m ca 0 
Mh ieee (76) 
3m 
0 0 i6 





B. STRUCTURAL MASS MATRIX 


For a shell having mass uw per unit circumferential length 


the kinetic energy in an axial length dz is taken to be 
t= sides woot vy) aaa . (77) 
0 


But, from equations (58) and (59), 


u, cos @ + Uy cos 26 ‘ (78) 


ce 
1) 


. > at ° 
and y= su, sin 6 =.5 u, sin 20 a (79) 
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where, from equations (62) and (63), 


e es e e TZ 

uy = =a) =o, sin = : (80) 
ee 2 1Z 

Uy = -43 sin mas (81) 


Making the appropriate substitutions and evaluating the 


resulting integrals over 6 and the length gives 


Soe 


4 
pau SMe 


FS) 
len 
Me 
ine) 
ine) 
+ 
Cu 
iS 
Oe 
Lo) 
NO 
v 
“- 
[ee] 
ine) 
Ww 


eT. =m m. 9445 + 
where Mee eTUAr 


= shell mass (including stiffening rings). 


Adding the bulkhead mass m_ to the rigid body mode and 


b 
including the added-mass-matrix the total kinetic energy 


becomes 


es e 2) 4 e e 
2T = (m + mM. + my, ) qy + = (m + m., ) 9495 


(83) 
1 are al aan 
+ 5} (m + m.) d5 + 6h (12m + 15m, ) 43 
Mm. 
For neutral buoyancy: m= m, +m, = het. B= ee. then 
equation (83) can be written as 
. 7 . 
ems take boiddiqd a (84) 
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where 


2 
2 aes 2, 0 
[M] = m |-£ (148) = > (148) 0 (85) 
0 0 ay (124158) 


C. STRUCTURAL STIFFNESS MATRIX 
The coefficients of the structural stiffness matrix are 
determined from energy considerations. The strain energy V 


is expressed in the form 
“4 Sy 
Pits wikliq} . (86) 
Since there is no strain energy associated with the 


rigid body mode and the other modes are elastically uncoupled, 


the stiffness matrix takes the form 


0 0 0°] 
ie a 0 Kyo 0 2 (87) 
0 0 ee 


an Shear Beam Mode 
The modal displacements are given by equations (59). 
There is no ring bending and shell bending contributions 
are negligible, so that only membrane energy need be 


Ponganered. 1.6.4, €, = 0, ee. = 0 and 






<ESTAM 22S yibetieg. Usel 


2 a 


xityam ehent ls. lexudcyctes =e So asgabediteon Me 
4 ‘ a. 


7 V aise ntesdes. off . an tarehtedon yaderts wot: feat 


ern”) aH rt bees ore 









Al Wid (Yaw Covseisoess vanes Hints oh Sh ster oe 


‘ di P ‘ A i 
i , a ‘ ; 
i abamoar: Vileaéivente sen a: et acto iy 
= ; : ; 
—% r 
— y , be Ores a8 wes 


_ 


ie 
> ay >: 
Lae 
? 
ud 
ae! 
¢ 
Pu 
' <4 
Py th: j u ; ¢ 





>a z q, sin © cos eo 7 (88) 


The strain energy Via associated with the shear beam mode is 


given by 
2 27 5 
Me ft) GY. a ceede (89) 
0 O 
or 
eaten)" 
LU ahem Sipe. (90) 
2. Buckling Mode 
The strain energy of the buckling mode is determined 


in two parts: the shell bending energy ie. and the membrane 
energy Vin" The modal displacements are given by equations (60). 


The shell bending energy is given by 


av, = : : (D ba + 2vD kK Kg t+ Deka + 4Digk gn) hoa dédz , 
(91) 
where D, and De are the flexural rigidities given as 
a er iseee oy, 
and 
Dy = a De) (93) 


af 


a 


(a9) , ae ae. 'S nik gp 7. -_ 


ai aben megad waade of) AdkW bat Bineeas 


168) , shee a ut 724 


(90) 
- 
frit Pata ah mF * 5 y . ye a0, 
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The value ® is a constant of proportionality determined by 
considering the orthotropic properties of the ring reinforced 


Saructure. Whe quantity D is the 'twisting' rigidity 


OZ 
given by 


D = D 3 . (94) 





where effects of ring twisting are neglected. 


The curvatures Ki» Ko and Kg are determined from equations 


(60) as: 


2 2 
ao ane : 5 ? Sele 5 [ 28 tu ? 
OZ a 36 
ae (95) 
és ewe Ih 3-u 
6z a d60z 





Substituting and performing the indicated integrations yields 


TD 3 
rates, Aare te a 2a, 3vTa 2 
eVh = ae [et Re ote "5 a = + BCL v)T z + Z 7 I3 
(96) 
The membrane energy is 
Lae T 5 5 
= 1 
ay. : : pete + GY 9 ) hga dédz : (97) 
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where 


Bee se0 eg Soh yo eos 28, (98) 
= Oy, Lew . cob ee 2m 
aes q3 ae 8! eee ee (99) 


Performing the indicated substitutions and evaluating the 


integrals gives 


math. Nee 








2 SUS a Boe ek ye 2 
eV, 5 Fi mY EEK aE + 5 64 3 (100) 
The parameter A is determined by minimizing the 
oV 
oO 
membrane energy with respect to A, i.e., a =O) anne 
result of this minimization is 
Dep ye 
is 
A = 5 a 4 (101) 
! 
h(n? i + 43) 
a 


The required coefficients of the stiffness matrix [K] 


are determined from equations (90), (96), (100) and (101). 


D. GENERALIZED FORCES 
The generalized force vector {f} is determined by employing 
the virtual work principle. Let 6W be the work done by the 


pressure p during a virtual displacement 6q. Then 


Bg 


au du au 
6bW=-/S ff p (— 6q, + —— 6q, + —— Sq.) a dzde@ (102) 
0 0 Leslie 2 ag aaa 
where au . gu and es are determined from the modal 
8a, * 345 243 
displacement equations (58-60), i.e., - = = cos 6, 
© 
oo sae sin ™ cos 6 5 oe 2G = sin® 22 cos 26 |. Evaluating 
845 R 3 Q 
the integral over the surface 
6W = [8q,,645,64,] [f,] (103) 
‘2 
9 
3 
where 
oT 
a ek J) p Gon 8} de, (104) 
0 
oT 
i, = eat ij pcos 6dée , (105) 
0 
al oT 
f,=-=—/s p cos 26 dé , (106) 
3 2 0 


The generalized forces are calculated from the two- 
dimensional finite element pressure-time history. The 
resulting generalized forces are given in Fig. 14. The 


rigid body force f, is in close agreement with the results 


ik 
calculated by Sette [1] and Murray [3] for a step pulse. 
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fy - rigid body force 


f, ~- shear beam force 


fr. - buckling force 





ct/a 


Fig. 14. Generalized Forces-vs.-Time 
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E. GENERALIZED DISPLACEMENTS 

The Newmark-8 method, equation (42), of time integration 
was employed to evaluate eee time dependent response of the 
structure. Figure 15 shows the resulting generalized 
displacements (q,E/p*a) as a function of the nondimensional 


time parameter ct/a. 


F, STRUCTURAL STRESSES 

The structural stresses are determined from the strain 
relations utilized to determine the coefficients of the 
Structural stiffness matrix. 


For the shear beam mode equation (88) gives 


=Gy,=G7 


a Z 43 sin ® cos 3 . (107) 


For the buckling mode the stress due to ring bending is 





2 6 ac are 
3Ec 
= 5 43 sin® -= eos 20 7; (108) 


where Cy is the distance from the centroid of the cross- 


section to the extreme fiber. The membrane tension stress 


(equation (98)) is 


Ss ¢@ = hte. = opr TA dz cos ee ese (109) 
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and the membrane stress (equation (99) is 


where X is determined from equation (101). 


Table I shows the maximum nondimensional stresses S,/p* 


as a funetion of ct/a. 


MAXIMUM MODAL STRESSES 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The results of this investigation confirm that the finite 
element method provides a satisfactory numerical method of 
evaluating the interaction of a plane acoustic blast wave 
and a rigid structure. Extension of the method to other two 
and three-dimensional geometries using the techniques 
employed appears to be practical. | 

The submarine structural model is a simple three degree- 
of-freedom model. It was chosen to provide an illustrative 
example of the implementation of the superposition theorem 
and not as a general model of a submarine. 

It is recognized that a more sophisticated model is 
necessary to determine the response of a submarine subjected 
to blast wave loading. To improve the model it would be 
advantageous to: 


a. Include additional modes, particularly a 'breathing' 
mode to incorporate the effects of hydrostatic loading; 


b. Distribute and elastically mount the mass of the sub- 
marine not included in the shell; 


c. Remove the restrictions on the rigidity of the 
bulkheads. 


Incorporating the above improvements and nonlinear 
structural characteristic into the submarine model will 
provide for a more complete analysis of the dynamic response 


of a submarine to blast-wave loading. 
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APPENDIX A 


PRESSURE-vs.-TIME AT THE RIGID CYLINDER-FLUID 
INTERFACE FOR 6 = 30°, 60°, 90°, 120° and 150° 


The pressure-time histories at the rigid cylinder-fluid 
interface at 6 = 30°, 60°, 90°, 120° and 150° are shown in 
Figs. 16-20. The solid line represents Sette's [1] results 
Gisplaced along the time axis to coincide with the incidence 
of the constant portion of the cubic ramp-step wave. 

Sette [1] provided no data for ® = 120° and ® = 150° and 
therefore no comparison is available. Data for ® = 0° and 
6 = 180° is presented in Figs. 9 and 10. 


The integration, wave and mesh characteristics are: 


Segment 
a 6 > 
et 
ra ae 
pal 
a Pie 
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Fig. 16. Pressure-vs.-Time 
Fluid Interface ( 
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Fig. 17 Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (6 = 60°) 








seb edna vanes er——an iol 





te 





ct/a 


Fig. 18. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (6 = 90°) 
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Fig. 19. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (8 = 120°) 
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Fig. 20. Pressure-vs.-Time at the Rigid Cylinder- 
Fluid Interface (8 = 150°) 
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APPENDIX B 


TABULATED PRESSURE-TIME HISTORY 
AT THE RIGID CYLINDER-FLUID INTERFACE 


The following tabulated data is the result of a two- 
dimensional finite element solution of the interaction of 
a ramp-step wave and a rigid cylindrical structure. 


The integration, wave and mesh characteristics are: 


h = a/6c, 

ets 0.97 a/c, 
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APPENDIX C 


SUBMARINE STRUCTURAL PARAMETERS 


The specific structural parameters used in the calculation 


of the coefficients of equation (57) are: 


p = fluid density = 9.36 x 107? lbf sec*/in’ 
=a Ok 10° ps2’, 

GPa, Gx 10° pel, 

v = 0.300, 

a '=)\hull radius = 180 tn., 

2 = bulkhead spacing = 480 in., 

On ies shell thickness = 1.60 in., 

B = m ,/m = 0.200 


The internal ring reinforcing frames are characterized 


by: 


i 


Frame spacing |e Ver Ula & oe 
Web thickness = 1.10 in., 
Frame depth = 9.20 in., 
Flange width = 9.20 in., 
Flange thickness = 1.80 in., 


ao = orthotropic constant of proportionality = 127. 
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The mass matrix is 


9.405 3.592 0. 
[M] = | 3.592 2,821 0. % 10° IWF sec*/in. 
0. 0. 1.102 


The stiffness matrix is 


0 0. 0 
[K] = 0. 107.9 0. x 10° 1bf/in. 


0. 0. 213236 
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